library(tidyverse)
library(ggsci)
library(ggthemes)
library(Cairo)
library(sf)
library(grid)
library(gridExtra)
library(ggrepel)
library(scales)
library(broom)
library(margins)
library(Hmisc)

#Abre base ---------------------------
rm(list=ls())
latinobarometro <- haven::read_dta("latinobarometro20062020.dta")

latinobarometro$pais_f <- factor(latinobarometro$idenpa, labels = 
                                             c("Arg", "Bol", "Bra", "Chi", "Col", "Cri", "Dom", "Ecu", 
                                               "Slv", "Gua", "Hon", "Mex", "Nic", "Pan", "Par", "Per", 
                                               "Uru"))

latinobarometro$pais_largo_f <- factor(latinobarometro$idenpa, labels = 
                                   c("Argentina", "Bolivia", "Brasil", "Chile", "Colombia", "Costa Rica", 
                                     "Rep. Dominicana", "Ecuador", 
                                     "El Salvador", "Guatemala", "Honduras", "México", "Nicaragua", 
                                     "Panamá", "Paraguay", "Perú", "Uruguay"))

latinobarometro$sex_f <- as_factor(latinobarometro$sex, levels = "labels")

#Gráficos-----------

tabla_medias <- latinobarometro %>% 
  group_by(pais_largo_f, numinves) %>% 
  summarise(media = weighted.mean(subh, wt, na.rm = TRUE), sd = sqrt(wtd.var(subh, wt, na.rm = TRUE)), se = sd / sqrt(length(subh))) %>% 
  ungroup() %>% 
  add_row(pais_largo_f ="Argentina", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Bolivia", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Brasil", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Chile", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Colombia", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Costa Rica", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Rep. Dominicana", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Ecuador", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="El Salvador", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Guatemala", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Honduras", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="México", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Nicaragua", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Panamá", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Paraguay", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Perú", numinves = 2012, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Uruguay", numinves = 2012, media = NA, sd = NA) %>%
  
  add_row(pais_largo_f ="Argentina", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Bolivia", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Brasil", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Chile", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Colombia", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Costa Rica", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Rep. Dominicana", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Ecuador", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="El Salvador", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Guatemala", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Honduras", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="México", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Nicaragua", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Panamá", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Paraguay", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Perú", numinves = 2014, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Uruguay", numinves = 2014, media = NA, sd = NA) %>%
  
  add_row(pais_largo_f ="Argentina", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Bolivia", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Brasil", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Chile", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Colombia", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Costa Rica", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Rep. Dominicana", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Ecuador", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="El Salvador", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Guatemala", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Honduras", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="México", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Nicaragua", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Panamá", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Paraguay", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Perú", numinves = 2015, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Uruguay", numinves = 2015, media = NA, sd = NA) %>%
  
  add_row(pais_largo_f ="Argentina", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Bolivia", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Brasil", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Chile", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Colombia", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Costa Rica", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Rep. Dominicana", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Ecuador", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="El Salvador", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Guatemala", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Honduras", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="México", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Nicaragua", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Panamá", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Paraguay", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Perú", numinves = 2016, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Uruguay", numinves = 2016, media = NA, sd = NA) %>%
  
  add_row(pais_largo_f ="Argentina", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Bolivia", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Brasil", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Chile", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Colombia", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Costa Rica", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Rep. Dominicana", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Ecuador", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="El Salvador", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Guatemala", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Honduras", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="México", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Nicaragua", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Panamá", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Paraguay", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Perú", numinves = 2017, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Uruguay", numinves = 2017, media = NA, sd = NA) %>%
  
  add_row(pais_largo_f ="Argentina", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Bolivia", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Brasil", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Chile", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Colombia", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Costa Rica", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Rep. Dominicana", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Ecuador", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="El Salvador", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Guatemala", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Honduras", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="México", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Nicaragua", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Panamá", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Paraguay", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Perú", numinves = 2019, media = NA, sd = NA) %>% 
  add_row(pais_largo_f ="Uruguay", numinves = 2019, media = NA, sd = NA)

#tabla_promedio <- latinobarometro %>% 
 # group_by(numinves) %>% 
#  summarise(media = mean(subh, na.rm = TRUE), sd = sd(subh, na.rm = TRUE), se = sd / sqrt(length(subh))) %>% 
 # mutate(pais_f = "Promedio países")

#tabla_medias <- tabla_medias %>% 
 # bind_rows(tabla_promedio)

#rm(tabla_promedio)

ggplot(tabla_medias, aes(x= as.factor(numinves), y=media, group = 1)) +
  geom_hline(yintercept = 5, color = "red", size = 0.5, linetype = "dashed") +
  geom_line(size = 0.5) +
  geom_pointrange(aes(ymin = media-sd, ymax = media+sd), size = 0.3) +
  labs(y="Promedio del estatus subjetivo") +
  theme_dark() +
  theme(axis.text.x = element_text(angle = 90, vjust=0.5, size = 8),
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 9)) +
  scale_y_continuous(limits = c(0,10), breaks = seq(0, 10, 2)) +
  facet_wrap(~pais_largo_f)


ggsave(filename = "grafico_sub_pais_anio.png", dpi = 300, type = "cairo", width = 8, height = 6)



##Medias del estatus objetivo--------------
medias <- latinobarometro %>% 
  group_by(pais_f) %>% 
  summarise(media = mean(status_obj, na.rm = TRUE)) %>% 
  arrange(desc(media))

latinobarometro %>% 
  group_by(numinves) %>% 
  summarise(media = mean(subh, na.rm = TRUE)) %>% 
  arrange(desc(media))

cuadro <- latinobarometro %>% 
  group_by(numinves, pais_f) %>% 
  summarise(media = mean(subh, na.rm = TRUE)) %>%
  spread(pais_f, value = media)

cuadro_sd <- latinobarometro %>% 
  group_by(numinves, pais_f) %>% 
  summarise(sd = sd(subh, na.rm = TRUE)) %>% 
  arrange(desc(sd))


##Mapas---------

subjetiva <- latinobarometro %>% 
  mutate(posiciones_medias_sub = case_when(subh >= 4 & subh <= 6 ~ 1,
                                      subh < 4 | subh > 6 ~ 0)) %>% 
  group_by(pais_f) %>% 
  summarise(media_sub = round(weighted.mean(subh, wt, na.rm = TRUE), digits = 2),
            media_obj = round(weighted.mean(status_obj, wt, na.rm = TRUE), digits = 2),
            media_gini = round(weighted.mean(gini, wt, na.rm = TRUE), digits = 2),
            media_gini_wid = round(weighted.mean(gini_wid, wt, na.rm = TRUE), digits = 2),
            media_gdp = round(weighted.mean(log_gdp, wt, na.rm = TRUE), digits = 2),
            media_gdp_wid = round(weighted.mean(log_gdp_wid, wt, na.rm = TRUE), digits = 2),
            media_sit = round(weighted.mean(situ_econ, wt, na.rm = TRUE), digits = 2),
            media_p10 = round(weighted.mean(p10, wt, na.rm = TRUE), digits = 2))


america <- st_read("América.shp")
america <- america %>% 
  filter(!PAÍS %in% c("Canadá", "Groenlandia (Dinamarca)", "Estados Unidos", 
                      "San Pedro y Miquelón (Francia)"))

merge_america <- merge(america, subjetiva, by="pais_f", all.x = TRUE)

mapa_sub <- ggplot() +
  geom_sf(data = merge_america, aes(fill = media_sub), color = "black", size = .2) +
  theme_map() +
  labs(fill = "Promedio estatus \nsubjetivo") +
  theme(legend.title = element_text(size = 8, face = "bold"),
        legend.text = element_text(size = 7),
        legend.key.size = unit(.85,"line")) +
  scale_fill_viridis_c(option = "inferno", direction = -1, na.value = NA, guide = "colourbar",
                       limits = c(3,5), breaks = seq(3, 5, 0.5)) 

mapa_obj <- ggplot() +
  geom_sf(data = merge_america, aes(fill = media_obj), color = "black", size = .2) +
  theme_map() +
  labs(fill = "Promedio estatus \nobjetivo") +
  theme(legend.title = element_text(size = 8, face = "bold"),
        legend.text = element_text(size = 7),
        legend.key.size = unit(.85,"line")) +
  scale_fill_viridis_c(option = "inferno", direction = -1, na.value = NA, guide = "colourbar",
                       limits = c(2.5,5), breaks = seq(2.5,5, .5)) 

#mapa <- grid.arrange(mapa_sub, mapa_obj, ncol = 2,
 #                    bottom = textGrob("Fuente: Elaboración propia en base a Latinobarómetro 2006-2013",
  #                   gp = gpar(fontsize = 9),
   #                  hjust = 1,
    #                 x = 1))

mapa <- grid.arrange(mapa_sub, mapa_obj, ncol = 2)

ggsave(mapa, file = "mapas.png", width = 6, height = 4.5, type = "cairo-png")




##Gráficos de medias------------
ggplot(subjetiva, aes(x= media_obj, y=media_sub, label=pais_f)) +
  geom_smooth(linetype = "twodash", color = "red", size = .5, method="lm", se=T) +
  geom_point(size = .7) +
  geom_text_repel(size = 3) +
  labs(y="Estatus subjetivo",
       x="Estatus objetivo") +
  theme_dark() +
  theme(title = element_text(size = 9),
        axis.text.x = element_text(size = 9),
        axis.text.y = element_text(size = 9),
        axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 10)) +
  scale_y_continuous(breaks = pretty_breaks(n = 8)) +
  scale_x_continuous(breaks = pretty_breaks(n = 8))

ggsave(file = "subj_obj.png", width = 6, height = 4.5, type = "cairo-png")


media_pib <- ggplot(subjetiva, aes(x= media_gdp_wid, y=media_sub, label=pais_f)) +
  geom_smooth(linetype = "twodash", color = "red", size = .5, method="lm", se=T) +
  geom_point(size = .7) +
  geom_text_repel(size = 2) +
  labs(title = "A",
       x="PIB per cápita (log)",
       y="Estatus subjetivo") +
  theme_dark() +
  theme(title = element_text(size = 9),
        axis.text.x = element_text(size = 7),
        axis.text.y = element_text(size = 7),
        axis.title.x = element_text(size = 7),
        axis.title.y = element_text(size = 7)) +
  scale_y_continuous(breaks = pretty_breaks(n = 8)) +
  scale_x_continuous(breaks = pretty_breaks(n = 8))


media_gini <- ggplot(subjetiva, aes(x= media_gini, y=media_sub, label=pais_f)) +
  geom_smooth(linetype = "twodash", color = "red", size = .5, method="lm", se=T) +
  geom_point(size = .7) +
  geom_text_repel(size = 2) +
  labs(title = "B",
       y="Estatus subjetivo",
       x="Gini") +
  theme_dark() +
  theme(title = element_text(size = 9), 
        axis.text.x = element_text(size = 7),
        axis.text.y = element_text(size = 7),
        axis.title.x = element_text(size = 7),
        axis.title.y = element_text(size = 7)) +
  scale_y_continuous(breaks = pretty_breaks(n = 8)) +
  scale_x_continuous(breaks = pretty_breaks(n = 8))

media_p10 <- ggplot(subjetiva, aes(x= media_p10, y=media_sub, label=pais_f)) +
  geom_smooth(linetype = "twodash", color = "red", size = .5, method="lm", se=T) +
  geom_point(size = .7) +
  geom_text_repel(size = 2) +
  labs(title = "C",
       y="Estatus subjetivo",
       x="Participación en el ingreso del top 10%") +
  theme_dark() +
  theme(title = element_text(size = 9),
        axis.text.x = element_text(size = 7),
        axis.text.y = element_text(size = 7),
        axis.title.x = element_text(size = 7),
        axis.title.y = element_text(size = 7)) +
  scale_y_continuous(breaks = pretty_breaks(n = 8)) +
  scale_x_continuous(breaks = pretty_breaks(n = 8))


media_sit <- ggplot(subjetiva, aes(x= media_sit, y=media_sub, label=pais_f)) +
  geom_smooth(linetype = "twodash", color = "red", size = .5, method="lm", se=T) +
  geom_point(size = .7) +
  geom_text_repel(size = 2) +
  labs(title = "D", 
       x="Percepción negativa de la situación económica del país (en %)",
       y="Estatus subjetivo") +
  theme_dark() +
  theme(title = element_text(size = 9), 
        axis.text.x = element_text(size = 7),
        axis.text.y = element_text(size = 7),
        axis.title.x = element_text(size = 6),
        axis.title.y = element_text(size = 7)) +
  scale_y_continuous(breaks = pretty_breaks(n = 8)) +
  scale_x_continuous(breaks = pretty_breaks(n = 8))


medias_group <- grid.arrange(media_pib, media_gini, media_p10, media_sit, ncol = 2)



ggsave(medias_group, file = "medias_grupos.png", width = 6, height = 4.5, type = "cairo-png")


##Gráficos regresión -----
latinobarometro$pais_anio <- paste(latinobarometro$pais_f, latinobarometro$numinves, sep = " ")

coeficientes <- latinobarometro %>% 
  group_by(pais_anio) %>% 
  do(tidy(lm(subh ~ status_obj, data = .))) %>% 
  select(pais_anio, term, estimate) %>% 
  pivot_wider(names_from = term, values_from = estimate) %>% 
  rename(constante = "(Intercept)", pendiente = status_obj)

ggplot(coeficientes, aes(x= pendiente, y=constante, label = pais_anio)) +
  geom_point(size = 1.5) +
  geom_smooth(linetype = "twodash", color = "red", size = .8, method="lm", se=T) +
  geom_text_repel(data=subset(coeficientes, (pendiente < .15 & constante >= 4) | 
                                               (pendiente > .35 & constante < 2.5)), size = 3.2) +
  labs(y="Constante",
       x="Pendiente") +
  theme_dark() +
  theme(axis.text.x = element_text(size = 9),
        axis.text.y = element_text(size = 9),
        axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 10)) +
  scale_y_continuous(breaks = pretty_breaks(n = 8)) +
  scale_x_continuous(breaks = pretty_breaks(n = 8))

ggsave(file = "constantes.png", width = 6, height = 4.5, type = "cairo-png")


latinobarometro %>% 
  filter(pais_f %in% c("Bra", "Col", "Gua", "Nic")) %>% 
  ggplot(aes(x = status_obj, y = subh, group = pais_f)) +
    geom_smooth(method = "lm", size = .8, aes(color = pais_f)) +
    labs(y="Estatus subjetivo",
       x="Estatus objetivo",
       color = "Países") +
    theme_dark() +
    theme(axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 8),
        axis.title.x = element_text(size = 9),
        axis.title.y = element_text(size = 9)) +
    scale_y_continuous(breaks = pretty_breaks(n = 8)) +
    scale_x_continuous(breaks = pretty_breaks(n = 8)) +
    scale_color_nejm(labels = c("Brasil", "Colombia", "Guatemala", "Nicaragua"))

ggsave(file = "pendientes_paises.png", width = 6, height = 4.5, type = "cairo-png")


##Margins ----------

interaccion <- lm(subh ~ status_obj * pais_f, data = latinobarometro, weights = wt)
summary(interaccion)

int <- summary(margins(interaccion, variables = "status_obj", 
                       at = list(pais_f = c("Arg", "Bol", "Bra", "Chi", "Col", "Cri", "Dom", "Ecu", 
                                            "Slv", "Gua", "Hon", "Mex", "Nic", "Pan", "Par", "Per", 
                                            "Uru"))))

# ggplot(int, aes(x= pais_f, y=AME)) +
#   geom_hline(yintercept = .25, linetype = "dashed", color = "red", size = .5) +
#   geom_pointrange(aes(ymin = lower, ymax = upper), size = .4) +
#   labs(y="Efectos marginales medios") +
#   theme_dark() +
#   theme(axis.text.x = element_text(vjust=0.5, size = 9),
#         panel.grid.minor = element_blank(),
#         axis.title.x = element_blank(),
#         axis.title.y = element_text(size = 10),
#         axis.text.y = element_text(size = 9)) +
#   scale_y_continuous(limits = c(.1,.4), breaks = seq(.1, .4, .05)) 
# 
# 
# ggsave(filename = "grafico_int1.png", dpi = 300, type = "cairo", width = 6, height = 4.5)
# 
# 
# ggplot(int, aes(x= fct_reorder(pais_f, AME), y=AME)) +
#   geom_hline(yintercept = .25, linetype = "dashed", color = "red", size = .5) +
#   geom_linerange(aes(ymin = lower, ymax = upper), size = .4) +
#   labs(y="Efectos marginales medios") +
#   theme_dark() +
#   theme(axis.text.x = element_text(vjust=0.5, size = 9),
#         panel.grid.minor = element_blank(),
#         axis.title.x = element_blank(),
#         axis.title.y = element_text(size = 10),
#         axis.text.y = element_text(size = 9)) +
#   scale_y_continuous(limits = c(.1,.4), breaks = seq(.1, .4, .05)) 
# 
# 
# ggsave(filename = "grafico_int2.png", dpi = 300, type = "cairo", width = 6, height = 4.5)


ggplot(int, aes(x= fct_reorder(pais_f, AME), y=AME)) +
  geom_hline(yintercept = .38, linetype = "dashed", color = "red", size = .5) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = .5, width = .3) +
  labs(y="Efectos marginales medios") +
  theme_dark() +
  theme(axis.text.x = element_text(vjust=0.5, size = 9),
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 10),
        axis.text.y = element_text(size = 9)) +
  scale_y_continuous(limits = c(.2,.55), breaks = seq(.2, .55, .05)) 


ggsave(filename = "grafico_int3.png", dpi = 300, type = "cairo", width = 6, height = 4.5)

### Modelos 

library(lme4)
library(lmtest)
library(lmerTest)
library(estimatr) # Para calcular ICs bootstrap

# Carregar dados (substitua "seu_arquivo.csv" pelo nome do seu arquivo de dados)
dados <- haven::read_dta("latinobarometro20062020.dta")

names(dados)

# Centrar variáveis
dados <- dados %>%
  group_by(countryyear) %>%
  mutate(
    mediastatus_obj = mean(status_obj),
    status_obj_cent = status_obj - mediastatus_obj,
    mediaage = mean(age),
    age_cent = age - mediaage,
    age_cent2 = age_cent * age_cent,
    media_ideologia = mean(ideologia),
    ideologia_cent = ideologia - media_ideologia,
    log_gpd_mean = mean(log_gdp),
    log_gdp_cent = log_gdp - log_gpd_mean,
    gini_mean = mean(gini),
    gini_cent = gini - gini_mean,
    p10_mean = mean(p10),
    p10_cent = p10 - p10_mean,
    situ_econ_mean = mean(situ_econ),
    situ_econ_cent = situ_econ - situ_econ_mean,
    status_pib = status_obj_cent * log_gdp_cent,
    status_gini = status_obj_cent * gini_cent,
    status_econ = status_obj_cent * situ_econ_cent,
    status_p10 = status_obj_cent * p10_cent
  )


# Modelos REML
m1 <- lmer(subh ~ 1 + (1 | country) + (1 | countryyear), data = dados)
m2 <- lmer(subh ~ status_obj_cent + sex + age_cent + age_cent2 + (1 | country) + (1 | countryyear), data = dados)
m3 <- lmer(subh ~ status_obj_cent + sex + age_cent + age_cent2 + log_gdp_cent + gini_cent + p10_cent + situ_econ_cent + (1 | country) + (1 | countryyear), data = dados)
m4 <- lmer(subh ~ status_obj_cent + sex + age_cent + age_cent2 + log_gdp_cent + gini_cent + p10_cent + situ_econ_cent + status_pib + status_gini + status_p10 + status_econ + (1 | country) + (1 | countryyear), data = dados)
m5 <- lmer(subh ~ status_obj_cent + sex + age_cent + age_cent2 + log_gdp_cent + gini_cent + p10_cent + situ_econ_cent + status_pib + status_gini + status_p10 + status_econ + (1 | country) + (1 | countryyear) + (status_obj_cent | countryyear), data = dados)

# Teste de razão de verossimilhança
lrtest(m4, m5)

# Estimar os modelos usando bootstrap para obter ICs
# Lembre-se de instalar o pacote estimatr se ainda não estiver instalado: install.packages("estimatr")
m1_boot <- lm_robust(subh ~ 1, clusters = country, data = dados)
m2_boot <- lm_robust(subh ~ status_obj_cent + sex + age_cent + age_cent2, clusters = country, data = dados)
m3_boot <- lm_robust(subh ~ status_obj_cent + sex + age_cent + age_cent2 + log_gdp_cent + gini_cent + p10_cent + situ_econ_cent, clusters = country, data = dados)
m4_boot <- lm_robust(subh ~ status_obj_cent + sex + age_cent + age_cent2 + log_gdp_cent + gini_cent + p10_cent + situ_econ_cent + status_pib + status_gini + status_p10 + status_econ, clusters = country, data = dados)
m5_boot <- lm_robust(subh ~ status_obj_cent + sex + age_cent + age_cent2 + log_gdp_cent + gini_cent + p10_cent + situ_econ_cent + status_pib + status_gini + status_p10 + status_econ + status_obj_cent, clusters = country, data = dados)

# Calcular os intervalos de confiança (ICs) bootstrap
confint(m1_boot)
confint(m2_boot)
confint(m3_boot)
confint(m4_boot)
confint(m5_boot)





